Method and apparatus for locating the source of an unknown signal

ABSTRACT

A method of determining the location of an unknown source  10  transmitting a signal to satellite relays  14  and  16  comprises receiving replicas of the signal from the relays at receivers  18.  The receivers  18  also transmit and receive reference signals via respective relays  14  and  16.  All signals are downconverted, digitized and correlated with one another in pairs using a correlation function including a term which compensates for time varying differential frequency offset (DFO). Compensation for time varying differential time offset or time dilation is achieved by replicating or adding to signal samples and applying phase corrections. This procedure enables a correlation maximum and associated measurement results to be obtained despite the effects of relay satellite motion which mitigate against this. Results are used in a prior art geometrical technique to locate the unknown transmitter.

[0001] This invention relates to a method and apparatus for locating thesource of an unknown signal received by a plurality of signal relays.

[0002] In IEEE Trans. on Aerospace and Electronic Systems, Vol. AES-18,No. 2, March 1982, P C Chestnut describes the basic technique oflocating the source of an unknown signal such as a ground-basedcommunications antenna; it involves determining the time difference ofarrival (TDOA) and/or frequency difference of arrival (FDOA) of twosignals from the source relayed to receivers by respective interceptplatforms in the form of relay satellites in geostationary orgeosynchronous orbits. The signals are relayed by the satellites alongtwo independent signal paths to a receiving station, ieground-satellite-ground paths. One satellite lies in the main beam orlobe of the source antenna radiation pattern and the other in asidelobe. Each satellite accepts a signal (uplink) from the source,frequency shifts it using a turnaround oscillator and returns itsfrequency-shifted equivalent (downlink) to a ground receiver. The twosignal path lengths are normally unequal, and this gives two signalarrival times at the receiver differing by the TDOA value. FDOA is dueto relay satellite motion relative to the Earth and to one another,which Doppler shifts both downlink signal frequencies. The positions ofthe two satellites and the receiving station are known, and the locus ofpoints of constant TDOA or FDOA is in each case a surface whichintercepts the Earth's surface to define a curve referred to as a lineof position (LOP). Two measurements of TDOA or FDOA at different times,or one of each at one or more times, provides two LOPs which intersectat the position of the source to be located.

[0003] TDOA is also referred to as differential time offset (DTO) andFDOA as differential frequency offset (DFO) or differential Dopplershift.

[0004] The technique of determining TDOA and FDOA from two receivedsignals is described in IEEE Trans. on Acoustics Speech and SignalProcessing, Vol. ASSP-29, No. 3, June 1981 by S Stein in a paperentitled “Algorithms for Correlation Function Processing”. It is alsodescribed in U.S. Pat. No. 5,008,679 relating to a transmitter locationsystem incorporating two relay satellites and using both TDOA and FDOAmeasurements. The technique involves deriving the degree of correlationbetween the signals by multiplying them together and integrating theirproduct. Trial relative time shifts and frequency offsets are introducedin sequence between the signals and their correlation is determined foreach. The time shift and frequency offset which maximise the correlationare taken to be the required TDOA and FDOA, subject to correction forsignal propagation delays in satellite transponders and frequency shiftsin satellites and in processing.

[0005] The degree of correlation is determined from what is referred toby Stein as the cross ambiguity function or CAF A(τ,ν) defined by:$\begin{matrix}{{A\left( {\tau,v} \right)} = {\overset{\tau}{\int\limits_{0}}{{z_{1}^{*}(t)}{z_{2}\left( {t + \tau} \right)}e^{{- 2}\pi \quad i\quad v\quad t}{t}}}} & (1)\end{matrix}$

[0006] A(τ,ν) is the integral of the product of two signals z₁(t) andz₂(t) [complex or analytic versions of real signals s₁(t) and s₂(t)]after a trial time shift τ and a trial frequency shift ν have beenintroduced between them in processing after reception at the receivingstation. The asterisk in z₁*(t) indicates a complex conjugate. When amaximum value of the modulus of A(τ,ν), ie |A(τ,ν)|, is obtained, thisbeing a correlation peak in the surface |A(τ,ν)| as a function of thetwo variables τ and ν, the values of τ and ν for the peak are therequired TDOA and FDOA.

[0007] The system of U.S. Pat. No. 5,008,679 requires satellitepositions and velocities to be known accurately, and needs highly stablephase in ground station and satellite oscillators. It has bandwidthlimitations for satellite orbits inclined to the Earth's equatorialplane, and needs two receivers which are on the same site and havecommon time and frequency references.

[0008] International Patent Application No PCT/GB95/02211 under thePatent Co-operation Treaty, published as WO 97/11383, relates to atransmitter location system employing a reference signal passing via thesame satellite relays as the unknown signal and processed in phasecoherence with it. The reference signal is used to remove a number ofsources of error and limitations to which earlier techniques aresubject, giving improved accuracy and capability for use under a widerrange of conditions. Despite this improvement it has been foundsurprisingly that from time to time it can be impossible to discern acorrelation peak in the CAF surface |A(τ,ν)|—all that can be seen isnoise.

[0009] A related but different technique for counteracting sources oferror using a broad band approach is disclosed in U.S. Pat. No.5,594,452 to Webber et al.

[0010] It is an object of this invention to provide an alternativemethod and apparatus for transmitter location.

[0011] The present invention provides a method of locating the source ofan unknown signal received by a plurality of signal relays, the methodincluding the steps of:

[0012] (a) arranging for a plurality of receivers to receive replicas ofthe unknown signal via respective signal relays; and

[0013] (b) subjecting the replicas to correlation processing;

[0014] characterised in that correlation processing in step (b) isperformed with a complex correlation function (CCF) at least partlycompensated for change in the replicas' Differential Frequency Offset(DFO) with time due to relay motion relative to the source andreceivers.

[0015] In a preferred embodiment of the invention, correlationprocessing in step (b) is also performed with data sets adapted in phaseand subject to data replication or removal to counteract time dilationarising from relay motion relative to the source and receivers.

[0016] As indicated above it has been found that the reason for failureto obtain a correlation peak using a prior art ambiguity function is dueto relay motion relative to the source and receivers. This is verysurprising for geostationary satellite relays in particular, becausetheir motion has hitherto been treated as constant over a measurementinterval, which in turn implies that it affects all measurements at asite equally. Satellite motion alone would not have been expected tomake a correlation peak obtainable on some occasions but not others.Despite this, in accordance with the invention it has been discoveredthat signal correlation is affected by velocity and accelerationcomponents of the relay satellite along the lines joining it to thesource and receiver, which results in the replicas' DFO and DifferentialTime Offset (DTO) being time dependent. DFO variation can be compensatedas indicated above by adapting the correlation function, and wherenecessary DTO variation can be compensated also by adapting datasamples.

[0017] Correlation processing in step (b) may include introducing atrial time offset between the replicas and evaluating their correlationand iterating this to obtain a correlation maximum and derive at leastone of the replicas' DTO and DFO; it may be carried out with a CCFcontaining an exponent of a function of time having a first term whichis a constant DFO value and a second term which is a product of time anda constant value for rate of change of DFO with time, ie a constantdifferential frequency rate offset (DFRO) value, and wherein step (b)also includes introducing a trial value corresponding to DFRO prior toevaluating correlation and iterating one type of trial value for each ofthe other type, repeating for more values of the other type anddetermining a DFRO value appropriate to compensate at least partly forchange in DFO with time.

[0018] The CCF may be expressed as A(τ₀, b₁, b₂) given by:$\begin{matrix}{{A\left( {\tau_{0},b_{1},b_{2}} \right)} = {\overset{\tau}{\int\limits_{0}}{{z_{1}^{*}(t)}{z\left( {t + \tau_{0}} \right)}e^{{- 2}\pi \quad {i{({b_{1} + {2b_{2}t}})}}t}{t}}}} & (2)\end{matrix}$

[0019] where z₁ and z₂ are data sets representing two signals beingcorrelated after traversing different paths, the asterisk denoting thecomplex conjugate of z₁, T is a time over which data sets are taken, τ₀is the signals' DTO, and b₁ and b₂ are constants in a linearapproximation to the variation of their DFO with time t as follows:

DFO≡ν=b ₁+2b ₂ t   (3)

[0020] where b₁ is an initial constant DFO value at initiation of datarecordal and 2b₂ is DFRO as defined above.

[0021] The constant DFO value b₁ may be determined by Fouriertransformation of the Equation (2) product ofz₁*(t)z₂(t+τ₀)exp(−2πib₂t²) to the frequency domain (ie ignoring the b₁term) after each iteration of Step (b), b₁ being subsequently determinedas the frequency at which the CCF maximum occurs. This reduces thecomputation required because it avoids stepping through trial values ofb₁ and re-evaluating the CCF for each. For a data set of N samples itreduces the number of computations from of order N² to Nlog_(e)N, asubstantial reduction for large N˜10⁶.

[0022] Correlation may be performed with a sample data set adapted inphase and number of samples, samples being replicated in or removed fromthe set according which data set is selected for adaptation (z₁ or z₂)and also to the sign of the rate of change with time of DTO associatedwith the data set. Data samples in one of two data sets to be correlatedmay adapted in phase by samples following a removed sample beingmultiplied by a phase factor or by samples including and following areplicated sample being multiplied by such a factor, the phase factorbeing either e^(2πifΔt) or e^(2πifΔt) according to which of the two datasets is adapted and to the sign of DTO change with time, and wherein fis a signal band centre frequency of the samples after downconversionfor sampling and Δt is the interval between samples.

[0023] Samples may alternatively be selected for removal and replicationon the basis of the following equation: $\begin{matrix}{\tau_{m} = {{{{- \frac{\lambda}{c}}b_{1}^{\prime}t_{m}} - {\frac{\lambda}{c}b_{2}t_{m}^{2}}} = {m\quad \Delta \quad t}}} & (4)\end{matrix}$

[0024] where t_(m) is the time from initiation of sampling to the mthsample selected for alteration, m is the selection number, Δt is thesampling interval, τ_(m) is the time dilation (expressed as a number ofsampling intervals) and λ is the wavelength at the centre frequency ofthe signal sampling band before any frequency downconversion at areceiving site. The parameter b′₁ is an estimate of DFO (b₁) with anysystematic error (eg frequency shift at a satellite relay) counteractedby deriving an estimated correction using a reference signal. Theparameter b′₁ may alternatively be derived by using trial values of itand calculating the modulus of the CCF until a CCF maximum is obtained.It also is possible to exploit b₂ in the above expression as anindependent parameter, thereby dissociating time dilation and timevariation of DFO. Equation (4) is solved for t_(m) to determine times atwhich time dilation results in data sets z₁ and z₂ being out of step bysuccessive increments of Δt. If the effects of time dilation are not toosevere, a linear approximation may be employed by deleting the t_(m) ²term in Equation (4).

[0025] In an alternative aspect, the present invention provides alocation device for locating the source of an unknown signal received bya plurality of signal relays, and including:

[0026] (a) a plurality of receivers for receiving replicas of theunknown signal via respective signal relays;

[0027] (b) a correlation processor for correlating the replicas andobtaining a correlation maximum indicating at least one of their DTO andDFO;

[0028] characterised in that the correlation processor is arranged toperform correlation with a complex correlation function (CCF) at leastpartly compensated for change in DFO with time due to relay motionrelative to the source and receivers.

[0029] In a preferred embodiment of the invention, the correlationprocessor is also arranged to perform correlation with data sets adaptedin phase and subject to data replication or removal to counteract timedilation arising from relay motion relative to the source and receivers.

[0030] The correlation processor may be arranged to employ a CCFcontaining an exponent of a function of time having a first term whichis a constant DFO value and a second term which is a product of time anda constant value for rate of change of DFO with time, ie a constantdifferential frequency rate offset (DFRO) value, and wherein thecorrelation processor is also arranged to introduce trial values of DFROand to evaluate correlation for pairs of trial values of DTO and DFROiteratively to obtain a correlation maximum.

[0031] The correlation processor may be arranged to introduce a trialtime offset between the replicas, evaluate their correlation and iteratethis procedure to obtain a correlation maximum. It may employ a sampledata set adapted in phase and number of samples, samples beingreplicated in or removed from the set according which data set isselected for adaptation (z₁or z₂) and also to the sign of the rate ofchange with time of DTO associated with the data set. It may be arrangedto adapt samples in one of two data sets to be correlated in phase bymultiplying samples following a removed sample by a phase factor or bymultiplying samples including and following a replicated sample by sucha factor, the phase factor being either e^(2πifΔt) or e^(−2πifΔt)according to which of the two data sets is adapted and to the sign ofDTO change with time, and wherein f is a signal band centre frequency ofthe samples after downconversion for sampling and Δt is the intervalbetween samples.

[0032] The correlation processor may select samples for removal andreplication on the basis of the correlation function's time dimensionspread over an interval in which the samples were obtained. If a CCFpeak is obtainable from data without time dilation compensation, thetime over which correlation extends may be estimated. A selected dataset may be contracted or expanded by the removal or insertion of samplesaccording to whether data set z₁ or z₂ is adjusted and also according towhether time dilation is positive or negative. Alternatively, thecorrelation processor may select samples for removal and replication onthe basis of Equation (4) above or a linear approximation thereto

[0033] In order that the invention might be more filly understood, anembodiment thereof will now be described, by way of example only, withreference to the accompanying drawings, in which:

[0034]FIG. 1 illustrates signal propagation between an Earth-basedtransmitter, satellite relays and Earth-based receivers;

[0035]FIG. 2 indicates effects to be compensated by a location system ofthe invention;

[0036]FIG. 3 shows schematically processing circuitry of a locationsystem of the invention;

[0037]FIG. 4 is a two dimensional plot of the modulus of the prior artcross ambiguity function |A(τ,ν)| as a function of the two variables τand ν obtained using data affected by satellite motion relative tosource and receivers;

[0038]FIG. 5 is equivalent to FIG. 4 in all respects except that inaccordance with the invention the complex correlation function |A(τ, b₁,b₂)| is plotted.

[0039]FIGS. 6 and 7 are graphs of complex correlation function signalpower against frequency for equivalent situations except that in theformer case there is no relay motion compensation whereas in the lattercase there is.

[0040] FIGS. 8 to 11 are graphs of simulated noise-free results for DTOand DFO plotted against one another and showing time and frequencydistortion, frequency compensation, time compensation and time andfrequency compensation respectively.

[0041] FIGS. 12 to 14 are equivalent to FIGS. 8, 9/10 and 11respectively with addition of the effects of noise.

[0042] Referring to FIG. 1, an unknown transmitter 10 located in theUnited States of America 11 generates a signal which causes interferencewith satellite communications. The transmitter 10 has a radiationintensity pattern 12 with a sidelobe 12A directed to a first satellite14 in a geostationary orbit, to which its signal propagates along afirst uplink path l₁ ^(u) (dotted line). It also has a radiation patternmain lobe 12B directed to a second satellite 16 in such an orbit; itssignal propagates along a second uplink path l₂ ^(u) to the satellite 16and produces interference with communications signals using the latter.The unknown signal frequency is determined by spectrum analysisequipment which routinely monitors satellite channels. A typicalcommunications satellite operating at Ku band (11-14 GHz) has 16channels each 36 MHz wide and each capable of carrying 100communications signals. The superscript “u” to path references l₁ ^(u)and l² ^(u) denotes uplink from the unknown transmitter 10.

[0043] The satellites 14 and 16 are monitored by respective receivers inthe form of ground station antennas 18A and 18B located in England andreferred to collectively as receivers 18. They receive the signal fromthe unknown transmitter 10 and retransmit it along first and seconddownlink paths l₁ ^(d) and l₂ ^(d) to respective receivers 18. Here thesuperscript “d” denotes a downlink path to a monitoring receiver. Thereceivers 18 receive respective replicas of the unknown signal whichinclude path delays and Doppler frequency shifts from satellite motion.The receivers 18 also transmit reference signals to the satellites 14and 16 respectively, these signals being replicas of one another whichare phase-locked together when transmitted initially. The satellites 14and 16 return frequency shifted equivalents of the reference signals torespective receivers 18.

[0044] Referring now also to FIG. 2, in which elements previouslydescribed are like-referenced, effects to be counteracted by theinvention are indicated upon a schematic drawing ofsource/satellite/receiver geometry. These effects were discovered inexperiments to investigate the feasibility of processing large sets ofdata obtained using satellites in geosynchronous orbits. It was foundthat failures to obtain correlation peaks with acceptable signal tonoise ratio were due to specific aspects of satellite motion. This isvery surprising, because geostationary or geosynchronous satellitemotion relative to a source or receiver has hitherto been treated asbeing substantially constant over the relatively short sampling periodsrequired for geolocation of a source, which implies that it would havean equal effect on all measurements in that interval at an individualsite. One would not therefore expect satellite motion alone to beresponsible for making a correlation peak obtainable on some occasionsbut not others. Despite this, in accordance with the invention it hasbeen discovered that signal correlation is affected by satellitevelocity and acceleration components along the lines joining thesatellite to the source and receiver. As will be described later in moredetail, it has been found that during a geolocation measurement thesecomponents can be sufficiently different at different times to makecorrelation obtainable on some occasions but not others.

[0045]FIG. 2 is marked with arrows 20 to indicate parameters comprisingcomponents of satellite velocity and acceleration (eg v_(1u), a_(1u)) inthe directions of the source 10 and receivers 18, together with signaltime delay and associated signal paths as defined in the table below.The present invention involves processing which counteracts time varyingDFO and preferably also time varying time dilation; the latter isdepicted in FIG. 2 by a linear source time axis 22 and a non-linearreceiver time axis 24. The axes 22 and 24 indicate that there is a timevarying time separation or dilation between successive signals arrivingat a receiver 18 because of satellite motion. Time Delay τ_(xy), Path(τ_(xy)) or Direction Velocity Component (v_(xy), or a_(xy),) v_(xy) orAcceleration and Satellite Component a_(xy) or Receiver τ_(1u) Delay onUplink path 1₁ ^(u) to Satellite 14 τ_(1d) Delay on Downlink path 1₁^(d) to Receiver 18A. τ_(2u) Delay on Uplink path 1₂ ^(u) to Satellite16. τ_(2d) Delay on Downlink path 1₂ ^(d) to Receiver 18B. v_(1u)Satellite 14 velocity component along uplink path 1₁ ^(u) to Source 10.v_(1d) Satellite 14 velocity component along downlink path 1₁ ^(d) toReceiver 18A. v_(2u) Satellite 16 velocity component along uplink path1₂ ^(u) to Source 10. v_(2d) Satellite 16 velocity component alongdownlink path 1₂ ^(d) to Receiver 18B a_(1u) Satellite 14 accelerationcomponent along uplink path 1₁ ^(u) to Source 10. a_(1d) Satellite 14acceleration component along downlink path 1₁ ^(d) to Receiver 18A.a_(2u) Satellite 16 acceleration component along uplink path 1₂ ^(u) toSource 10. a_(2d) Satellite 16 acceleration component along downlinkpath 1₂ ^(d) to Receiver 18B.

[0046]FIG. 2 shows what would happen to points originally regularlyspaced on a waveform which passes over a path in space and is subjectedto time varying Doppler effect. An analogue to digital converter (ADC)is used to sample a received signal and gives a regular set of samplesto a signal which has undergone Doppler effect; these samples requiretiming adjustment, in one of the data sets z₁ or z₂ in Equation (2), inorder to correspond to samples in the other data set which haveundergone a different degree of Doppler effect. FIG. 2 does not portraythis—but instead what is shown is the result of differential Dopplereffect on an array of points on a waveform which were originallyregularly spaced. The dilation differs across four receiver channels tobe described later, and is shown schematically at 22 for a singlechannel.

[0047] Referring now also to FIG. 3, a signal processing system for atransmitter location system (TLS) of the invention is shownschematically and indicated generally by 30. It is emphasised that thisdrawing is schematic and does not show all system features in detail. Itwill be described in outline only, because transmitter location isdescribed in great detail in the prior art of U.S. Pat. No. 5,008,679and published PCT Appln No WO 97/11383. Readers of this specificationare referred to PCT Appln No WO 97/11383 for a detailed description ofhow transmitter location may be put into practice, except as regardsmeasurement of additional parameters and complex correlation functionprocessing to achieve satellite motion compensation in accordance withthis invention.

[0048] To simplify description, the processing system 30 is illustratedfunctionally rather than in terms of actual items of apparatus; someitems subsume multiple functions and others shown separately may beimplemented together as a single device. The processing system 30 hasfour channels 32AR, 32AU, 32BR and 32BU; suffixes A and B indicateassociation with receivers 18A and 18B respectively and suffixes R and Uindicate association with reference and unknown signals respectively.Components have like numerical references with differing suffixesdistinguishing channels. The channels and their components will bereferred to without suffixes to indicate any or all withoutdifferentiation, and with suffixes if necessary to be specific. Receiver18A is connected to channels 32AU and 32AR, and receiver 18B isconnected to channels 32BU and 32BR.

[0049] Each channel 32 contains a series-connected arrangement of amixer stage 34 for frequency down conversion and analogue filtering, ananalogue to digital converter (ADC) 36 and an exponentiation stage 38(H) for conversion of real signals to complex equivalents with in-phaseand quadrature components. The mixers 34 are connected to a common localoscillator signal feed LO. Outputs from exponentiation stages 38AU and38BU pass to a first correlation processing stage 40U, and that from theexponentiation stages 38AR and 38BR to a second correlation processingstage 40R. The processing stages 40U and 40R carry out parameterdetermination and correlation function processing as will be describedlater in more detail. Geometrical calculation and mapping to locate theunknown transmitter 10 is done in a mapping stage 42. Theexponentiation, correlation processing and mapping stages 38, 40 and 42may be implemented in a single computer.

[0050] The processing system 30 operates as follows. Signals from thereceivers 18 are amplified and fed to the mixer stages 34 where they arefiltered and downconverted to a centre frequency suitable for analogueto digital (A/D) conversion. Downconversion takes place in severalstages, with filtering to define the relevant signal frequency band ineach case and to prevent aliasing at the A/D stage. Mixer stages 34AUand 34BU are tuned to extract the unknown signal, and mixer stages 34ARand 34BR are tuned to extract the reference signal.

[0051] The ADCs 36 sample signals from the respective mixer stages 34using automatic level control (ALC) to minimise the effect ofquantisation noise. The signal samples so formed are real, and theexponentiation stages 38 convert them to analytic or complex form usingthe Hilbert Transform. In particular, a set s₁(t) of signal samples fromchannel 32AU is converted to an analytic sample set z₁(t) which is thenconjugated to form a set z₁*(t). A similar procedure but withoutconjugation is followed for channel 32BU to produce a sample set z₂(t).These two analytic sample sets for the unknown signal are thencorrelated together at 40U as will be described later to determine theirDFO and DTO. Like sample sets are produced at 38AR and 38BR for thereference signal and are correlated at 40R.

[0052] The real and imaginary parts of the complex sample are denoted asfollows: channel 32AU: (x₁(t_(i)), y₁(t_(i)) or z₁; channel 32BU(x₂(t_(i)), y₂(t_(i))) or z₂; channel 32AR (r₁(t_(i)), s₁(t_(i))) or p₁;channel 32BR (r₂(t_(i)), s₂(t_(i))) or p₂.

[0053] Each correlation processing stage 40 correlates two sets ofsignal samples obtained respectively from distorted versions of the sameoriginal signal, the distortion arising from propagation via satelliterelays; ie the two versions of the unknown signal in channels 32AU and32BU are correlated with one another, as are the two versions of thereference signal in channels 32AR and 32BR. Each of the sample sets alsohas a noise component which is statistically independent of the noisecomponent of sample sets from other channels.

[0054] As in the prior art, the procedure is to equalise the sample setsused in correlation, by introducing adjustable parameters into acorrelation operation until the correlation function is a maximum; theparameters are related to the source/satellite/receiver geometry. Thereference signal is used as in the prior art to remove noise which iscorrelated with that in the unknown signal, and the reference signalsource provides a reference point with respect to which the unknownsource location is determined. However, in the prior art eachcorrelation operation involved only fixed values of differentialfrequency offset (DFO or ν) and differential time offset (DTO or τ). Inthe present invention, an expression for time varying DFO is used tocompensate for satellite motion and if necessary time varying DTO iscounteracted. If satellite motion is not too severe it is possible toobtain adequate results by compensating for time varying DFO only. Thecorrelation operation is changed to implement this using a newcorrelation function different to that of Equation (1). The new functionis referred to as the complex correlation function or CCF to distinguishit from the earlier expression; it is suitable for counteracting timevarying DFO using geosynchronous satellites in orbits inclined to theecliptic. It does not counteract time varying DTO, which is a differentprocedure to be described later. The CCF is expressed as A(τ₀, b₁, b₂)given by Equation (2) above and repeated below for convenience:$\begin{matrix}{{A\left( {\tau_{0},b_{1},b_{2}} \right)} = {\overset{\tau}{\int\limits_{0}}{{z_{1}^{*}(t)}{z\left( {t + \tau_{0}} \right)}e^{{- 2}\pi \quad {i{({b_{1} + {2b_{2}t}})}}t}{t}}}} & (2)\end{matrix}$

[0055] where z₁ and z₂ are sets of data samples representing two signalsbeing correlated, the asterisk denotes the complex conjugate of z₁, T isthe time over which samples are taken, τ₀ is a time invariant componentof the differential time offset (DTO) between the signals because oftheir different paths, and b₁ and b₂ are constants in the Equation (3)expression for DFO:

DFO≡ν=b ₁+2b ₂ t   (3)

[0056] The constant b₁ is an initial value of DFO at the beginning ofdata recordal in a measurement; 2b₂ is the rate of change of DFO withtime, referred to as the differential frequency rate offset or DFRO.

[0057] Equation (3) expresses the discovery in accordance with theinvention that in the case of geosynchronous satellites the variation ofDFO over a measurement can often be assumed to be linear with time to agood approximation. If the variation of DFO with time is sufficientlynon-linear, higher powers of time t—ie b₃t², b₄t³ etc—may be included inthe expression for DFO in Equation (3), but these can be omitted undermost circumstances. They would be included if the motion of thesatellite intercept platform were to require it, but at the expense ofadditional processing. The parameters τ₀ and b₂ are determined byinserting trial values of each into Equation (2), evaluating the CCF anditerating this procedure to find the maximum value of the correlationsurface |A(τ₀,b₁,b₂)|, ie the modulus of the CCF as a function of τ₀,b₁, b₂. It is not necessary to iterate b₁, because this is determined byinspection from the position of the CCF peak in the frequency domain aswill be described later in more detail.

[0058] The correlation function A(τ₀,b₁,b₂) is applicable to situationswhere the satellite relays or intercept platforms are geostationary(GEO) and large data sets are required, leading to long data collectiontimes and changing DFO. It is also applicable in situations where othertypes of intercept platforms are used, eg aircraft, GEO-LEO (low Earthorbit) MEO-MEO (medium Earth orbit) and LEO-LEO satellite combinations,but the additional technique of time dilation compensation is thenrelevant as will be described later.

[0059] In more detail, the correlation procedure using Equation (2) isas follows. Values for the delay and DFRO parameters τ₀ and b₂ arechosen. Initially, these parameters are set to zero, and each issubsequently incremented through a predetermined range in discretesteps: each step is a multiple of a basic increment determined by thetime Δt between successive signal samples (that is, a multiple of theinverse of the sampling rate). The sampling rate is numerically equal totwice the sampling bandwidth set by low pass filters incorporated in themixers 34.

[0060] The basic increments in τ₀ and b₂ are denoted Δt and Δb₂respectively; Δt=1/2B_(s) and Δb₂=1/2T², where B_(s) and T arerespectively the sampling bandwidth and sampling period of the ADCs 36.

[0061] For each selected value of DTO τ₀, and DFRO b₂, a fast Fouriertransform (FFT) of the inner product of z₁*(t) z₂(t+τ₀) multiplied byexp(−2πib₂t²) is taken in order to convert it to the frequency domain.This enables the DFO b₁ to be determined automatically as that frequencyoffset at which the CCF function maximum occurs. All values of DFO areaccounted for in a single FFT operation, so it is unnecessary toincrement in steps Δb₁ and evaluate the CCF after each step. The FFT isrepeated for each set of values of (τ₀, b₂); ie first values of τ₀ andb₂ are selected and an FFT of z₁*(t)z₂(t+τ)exp(−2πib₂t²) is evaluated.This is repeated for a range of values of τ₀ in incremental steps ofΔt=1/2B_(s); b₂ is then incremented by one 0.001 Hz/second increment(this is a nominally convenient value) and the range of values of τ₀ andassociated FFT evaluations is repeated. The procedure is repeated untilall values of τ and b₂ in their respective ranges have been used. Foroperation in the SHF frequency band, from 7.9 GHz to 8.4 GHz, withGEOsat intercept platforms, a typical value of b₂ was found to be 0.05Hz/second. The modulus of the complex correlation function (CCF) has amaximum in the frequency domain, and the values of τ₀ and b₁ at thismaximum provide the required DTO and DFO. This provides acomputationally efficient method of evaluating the CCF for a range ofvalues of (τ₀, b₁, b₂).

[0062] The foregoing procedure locates the CCF maximum approximately.The location accuracy can be improved by interpolation as described indetail in International Patent Application No WO 97/11383 referred toabove in relation to the complex ambiguity function (CAF). A range ofvalues of the modulus (or modulus squared) of the CCF about its maximumvalue is taken, and interpolation is used to obtain a better estimate ofthe DTO, DFO and DFRO values.

[0063] This procedure results in values for τ₀, b₁ and b₂ for the signalfrom the source to be located. It is repeated for the reference signalto obtain corresponding values τ_(0,ref), b_(1,ref), and b_(2,ref). Theprocedure gives DFO corrected to the instantaneous value it had at thebeginning of data gathering when t=0; Equation (3) refers.

[0064] The values of τ₀ and τ_(0,ref) are used to determine the uplinktime difference of arrival TDOA_(u), ie the difference between the timesof arrival of the two signals or signal replicas from the unknown source10 at respective satellites 14 and 16. The values of b₁, b_(1,ref), b₂and b_(2,ref) are used to determine the DFO values for the unknownsignals and the reference signals. A theoretical treatment of this willnow be given.

[0065] Referring firstly to use of τ₀ and τ_(0,ref), DTO is thedifference in the time for two distorted versions of an original unknownor reference signal to traverse a two different paths to a receiver, ieeither from the source 10 to a receiver 18 or from and to a receiver 18via satellite 14 or 16. It can be represented as a function of time tover the duration of the sample, in terms of independent, time-constantparameters (τ₀, b₁, b₂, . . . etc) related to the kinematics of thesource/satellite/receiver system:

DTO≡τ=τ ₀ +a ₁ t+a ₂ t ² +a ₃ t ³+  (4)

[0066] The DTO τ is formed from the difference of two terms, τ₁ and τ₂,which represent signal propagation delays arising from traversal ofdifferent paths by replicas of an single original signal. It may beexpressed as the sum of uplink and downlink delays plus any additionaldelay experienced in traversing electronic signal channels (e.g. asatellite transponder system) as follows:

τ=τ₁−τ₂=(τ_(1u)+τ_(1d))−(τ_(2u)−τ_(2d))=(τ_(1u)−τ_(2u))+(τ_(1d)−τ_(2d))  (5)

[0067] where terms are defined in Table 1. Equation (5) represents theDTO in terms of the sum of the difference of uplink delays anddifference of downlink delays, and subsumes within these any differencebetween electronic delays.

[0068] Sampling periods required for geolocation are usually relativelyshort (milliseconds to tens of seconds), and despite satellite motion ithas been found that uplink and downlink delays can be represented astruncated power series in time t as follows:

τ_(1u)=τ_(0,1u) +a _(1,1u) t+a _(2,1u) t ²+  (6)

[0069] with similar expressions for other delays given in the tableearlier, ie τ_(1d), τ_(u2) and τ_(d2), with appropriate changes incoefficient indices to a_(1,1d) etc. The coefficients a_(1,1u) etc arenot the same as the accelerations defined in the earlier table, but canbe related to satellite velocities and accelerations. The fullexpression for the DTO becomes: $\begin{matrix}\begin{matrix}{\tau = \quad {\left( {\left( {\tau_{0,{1u}} - \tau_{0,{2u}}} \right) + \left( {\tau_{0,{1d}} - \tau_{0,{2d}}} \right)} \right) + \left( {\left( {a_{1,{1u}} - a_{1,{2u}}} \right) +} \right.}} \\{\quad \left( {a_{1,{1d}} - {\left. a_{1,{2d}} \right))t} + {\left( {\left( {a_{2,{1u}} - a_{2,{2u}}} \right) + \left( {a_{2,{1d}} - a_{2,{2d}}} \right)} \right)t^{2}} + \ldots} \right.}\end{matrix} & (7)\end{matrix}$

[0070] The parameters τ₀, a₁, and a₂ etc can be found by inspection, ieby equating coefficients of powers of time t in Equations (4) and (7).

[0071] Equation (7) for the DTO τ contains an uplink contribution whichrelates directly to the unknown source position to be determined, and adownlink contribution which is known from the satellite/receiver systemgeometry. The uplink contribution is referred to as the TDOA_(u). Inorder to remove biases from the estimate of TDOA_(u), an equivalentexpression for the DTO of the reference signal is subtracted from theDTO of the unknown source signal. The known geometry of thesource/satellite/receiver system allows the DTO of the reference signalto be calculated, apart from biases (e.g. satellite processing delays)affecting it in the same way as they affect the source signal DTO. Thisis why the reference signal DTO is also measured. The TDOA_(u) isrelated to the measured DTO values through the relation: $\begin{matrix}{{T\quad D\quad O\quad A_{u}} = {{\tau - \tau_{ref} + \frac{\left( {l_{1}^{u} - l_{2}^{u}} \right)_{ref}}{c}} = {\tau - \tau_{ref} + \frac{\left( {l_{1}^{d} - l_{2}^{d}} \right)_{ref}}{c}}}} & (8)\end{matrix}$

[0072] where c is the velocity of light and the subscript ‘ref’ in (8)indicates the paths are between the satellites and the receiving groundterminals over which the reference signal passes. DFO ν is thedifference in frequency between two received signals. It can beexpressed in terms of the kinematics of the source/satellite/receiversystem using the general relationship between Doppler, velocity and rateof change of time delay with time, exemplified for the satellite uplink,as follows: $\begin{matrix}{{\Delta \quad v_{1u}} = {{\frac{\quad}{t}\varphi_{1u}} = {{{- \frac{\quad}{t}}\left( \frac{l_{1}^{u}}{\lambda_{u}} \right)} = {\frac{v_{1u}}{\lambda_{u}} = {{- \frac{c}{\lambda_{u}}}\frac{\quad}{t}\tau_{1u}}}}}} & (9)\end{matrix}$

[0073] where Δν_(1u) is the Doppler shift in frequency occurring on theuplink path 1₁ ^(u) to satellite 14, τ_(1u) is the phase and λ_(u) thewavelength of a signal (at its centre frequency) as received at thesatellite 14 or 16 on uplink. Similar expressions for other systemuplinks and downlinks are obtainable by change of indices.

[0074] The unknown signal DFO is given by: $\begin{matrix}{{{DFO} \equiv v} = {{- {\frac{c}{\lambda_{u}}\left\lbrack {{\frac{\quad}{t}\tau_{1u}} - {\frac{\quad}{t}\tau_{2u}}} \right\rbrack}} - {\frac{c}{\lambda_{d}}\left\lbrack {{\frac{\quad}{t}\tau_{1d}} - {\frac{\quad}{t}\tau_{2d}}} \right\rbrack} + {bias}}} & (10)\end{matrix}$

[0075] where λ_(u) and λ_(d) are the free space wavelengths of theuplink and dowlink signals respectively at their centre frequencies, andbias arises from satellite turnaround frequency shift and from satellitepositional error.

[0076] The reference signal DFO is similarly given by: $\begin{matrix}{{{DFO}_{ref} \equiv v_{ref}} = {{- {\frac{c}{\lambda_{u,\quad {ref}}}\left\lbrack {{\frac{\quad}{t}\tau_{1u}} - {\frac{\quad}{t}\tau_{2u}}} \right\rbrack}_{ref}} - {\frac{c}{\lambda_{d,\quad {ref}}}\left\lbrack {{\frac{\quad}{t}\tau_{1d}} - {\frac{\quad}{t}\tau_{{2d} -}}} \right\rbrack}_{ref} + {bias}}} & (11)\end{matrix}$

[0077] Subtracting (11) from (10) cancels out bias which is the same inboth equations and gives an expression for the uplink frequencydifference of arrival FDOA_(u): $\begin{matrix}{{FDOA}_{u} = {v - v_{ref} + {\frac{c}{\lambda_{d}}\left\lbrack {{\frac{\quad}{t}\tau_{1d}} - {\frac{\quad}{t}\tau_{2d}}} \right\rbrack} - {\frac{c}{\lambda_{u,\quad {ref}}}\left\lbrack {{\frac{\quad}{t}\tau_{1u}} - {\frac{\quad}{t}\tau_{2u}}} \right\rbrack}_{ref} - {\frac{c}{\lambda_{d,\quad {ref}}}\left\lbrack {{\frac{\quad}{t}\tau_{1d}} - {\frac{\quad}{t}\tau_{2d}}} \right\rbrack}_{ref}}} & (12)\end{matrix}$

[0078] where the last three terms are known from the system geometry,and the first two terms are the unknown and reference DFOs measured bydetermining the CCF maximum as described earlier. A further measurableis also determined, the differential frequency rate offset (DFRO=2b₂),this being the difference between Doppler rates (DFO rates of change) ofsignals measured along respective paths from the source 10, through thesatellite intercept platforms 14 and 16, to the terrestrial interceptstations 18 as shown in FIG. 2. It may be used to form a parameter knownas the uplink differential Doppler rate of arrival (DDROA_(u)) whichexploits the acceleration effects of the source/satellite/receiversystem. The DDROA_(u) is the difference in Doppler-rates measured alongthe paths from the unknown source 10 to the satellite intercept platform14 and 16. It is useful for the location of moving sources. The DDROAcan be expressed in terms of the DFRO by direct differentiation of (12)as follows: $\begin{matrix}{{DDROA}_{u} = {{DFRO} - {DFRO}_{ref} + {\frac{c}{\lambda_{d}}\left\lbrack {{\frac{^{2}}{t^{2}}\tau_{1d}} - {\frac{^{2}}{t^{2}}\tau_{2d}}} \right\rbrack} - {\frac{c}{\lambda_{u,{ref}}}\left\lbrack {{\frac{^{2}}{t^{2}}\tau_{1u}} - {\frac{^{2}}{t^{2}}\tau_{2u}}} \right\rbrack}_{ref} - {\frac{c}{\lambda_{d,{ref}}}\left\lbrack {{\frac{^{2}}{t^{2}}\tau_{1d}} - {\frac{^{2}}{t^{2}}\tau_{2d}}} \right\rbrack}_{ref}}} & (13)\end{matrix}$

[0079] The uplink TDOA (TDOA_(u) as in (8)) and the uplink FDOA(FDOA_(u) as in (12)) are used in the process of emitter positionlocation of stationary RF sources (on the Earth's surface), as in theprior art of published International Application No WO 97/11383. In thisparticular case the DDROA_(u) is not used explicitly. However, even whenthis is the case, the effect of satellite acceleration, which gives riseto the DDROA_(u), is counteracted in the CCF as previously describedusing trial values of the parameter b₂ as described earlier. Thisparameter must be utilised in evaluating the CCF even when it is notneeded further in locating the source 10. This is because the RF emitterlocation problem has two parts. The first is detection of the peak ofthe CCF and the second is relating the peak of the CCF to a position onthe Earth's surface. The first phase is sensitive to the effects ofsatellite motion; requiring b₂ as a compensation parameter in additionto those used in prior art (τ₀ and b₁).

[0080] The source 10 is located by using the determined TDOA_(u) andFDOA_(u) values to the Earth's surface by prior art geometrical methods.Firstly, the parameters τ₀, b₁ and b₂ are obtained for the unknown andreference signals by maximising the CCF in Equation (2) as describedearlier. The values of b₁, b_(1,ref), b₂ and b_(2,ref) are used todetermine the signal and reference DFO values by insertion in Equation(3) with time t set to zero. This sets DFO as its instantaneous value atthe beginning of data gathering—it is convenient to correct DFO to thevalue at this time but not essential.

[0081] The values of ν and ν_(ref) are substituted into Equation (12) todetermine FDOA_(u). Here again time t is set to zero; a consequence ofthis is that the location of the source 10 is determined using satelliteephemeris data (instantaneous positions in space) associated with thetime at the start of sampling.

[0082] TDOA_(u) is multiplied by c the velocity of light to givedifferential slant range (DSR), this being the difference between thetwo signal path lengths to the signal relay satellites 14 and 16 from atransmitting station, ie either the reference transmitter location18A/18B or the unknown source 10. FDOA_(u) is multiplied by λ_(u) thesignal wavelength in the uplink path to give the differential slantrange rate (DSRR) or time rate of change of DSR. International PatentApplication No PCT/GB95/02211 under the Patent Co-operation Treaty,published as WO 97/11383, describes in detail how to locate the source10 on the Earth's surface from DSR and DSRR values using a Taylorexpansion approach, and will not be described here.

[0083] Counteracting only time varying DFO is sufficient for use withtwo geostationary satellite (GEOsat) relays when large sample sets arerequired. It focuses the CCF in the frequency dimension. Its effects canbe seen from FIGS. 4 and 5; FIG. 4 shows results obtained using theprior art ambiguity function of Equation (1) and FIG. 5 shows resultsobtained using the complex correlation function or CCF of Equation (2)in accordance with the invention. The same data sets were used in bothcases, obtained using two geostationary satellites. To meet patentrequirements, these drawings were produced by tracing from colour printsin which colour indicated magnitude, and because of this they are nothighly accurate. They do however provide a reliable indication of thebenefits to be obtained by use of the invention.

[0084]FIG. 4 is a pseudo-three-dimensional view of the prior art complexambiguity function of Equation (1) plotted against DTO and DFO. All thatcan be seen is noise. FIG. 5 is a similar plot of the complexcorrelation function of Equation (2) in accordance with the invention.In FIG. 5 a sharp correlation maximum 50 is fully resolved and standswell above the noise at 52.

[0085] Referring now to FIG. 6, the effect of satellite motion is shownonce more in a plot of CCF modulus squared or power against frequency.There should be a single CCF peak at the DFO of −820.38, but insteadthere is a plateau 60 spread out over more than 16 Hz. This iscontrasted with FIG. 7, in which DFO variation has been compensatedgiving a well resolved peak 70.

[0086] It may be necessary to use fast moving relay platforms forgeolocation, eg aircraft, medium Earth orbit satellites (MEOsat), lowEarth orbit satellites (LEOsat), or combinations of these withgeostationary satellites (GEOsat). If so the CCF must also be focused inthe time dimension to compensate for time dilation. This is necessary tocompensate for a time delay τ₀ that varies appreciably in the course ofdata collection for a measurement.

[0087] The CCF of Equation (2) only compensates for DFO variation in thesampling interval T. It does not compensate for time dilation, ie DTOvariation in this interval. DTO variation is counteracted by periodicreplication or excision of samples to compensate for samples in one dataset (z₁ or p₁) being out of step with those in the associated data set(z₂ or p₂). For a measurement using two satellites, one geostationaryand the other low earth orbit, it might involve inserting or removingone sample in every few hundred in a data set having in the order of amillion samples. A phase correction is also required, and is applied toall samples in a data set following an insertion or excision, togetherwith the inserted sample if this occurs. Successive phase correctionsaccumulate, each one being applied from the point of insertion orexcision through the subsequent samples in the data set.

[0088] The timing of the sample adjustments, insertions/excisions isdetermined from the time taken for a sample slip to occur. A slip occurswhen the value of time dilation changes by an amount equal to thesampling interval at the ADC. The time to the mth slip (total timedilation of m sampling intervals) is estimated by one of two approaches:either (a) from the spread of the ambiguity surface, in the timedirection, across the sampling interval, or (b) from use of the τ₀, b₁,b₂ parameters obtained from maximising the CCF as described earlier.

[0089] In approach (a), the correlation processor selects samples forremoval from or replication in a set of data samples on the basis of thecomplex correlation function's time dimension spread over an interval inwhich the data set was obtained. If a CCF peak can be obtained from datawithout time dilation compensation, the number of time intervals Δt overwhich the correlation is spread is estimated. For a spread of K timeintervals, the number of samples between positions of adjustment(successive sample slips) of the data set is N/K where N is the numberof samples in the data set. The data set z₁ or z₂ (see Equation (2))selected for adjustment is a matter of implementation convenience andtherefore arbitrary in principle; eg if z₁ corresponds to a signalexpanded by time dilation, it can be contracted by sample removal tocorrect it, or alternatively z₂ may be expanded to match it. Theselected data set is contracted or expanded by the removal or insertionof samples according to which of the two data sets is selected foradjustment and also according to whether the time dilation determined ispositive or negative; if z₁ is selected, it is expanded if the timedilation is negative, which corresponds to rate of change of DTO beingnegative.

[0090] In the alternative approach (b), the correlation processor 30selects samples for removal and replication on the basis of thefollowing equation, which is an approximation derived from satellitemotion and is quadratic in time: $\begin{matrix}{\tau_{m} = {{{{- \frac{\lambda}{c}}b_{1}^{\prime}t_{m}} - {\frac{\lambda}{c}b_{2}\tau_{m}^{2}}} = {{m\Delta}\quad t}}} & (14)\end{matrix}$

[0091] t_(m) is the time from initiation of sampling to the mthreplicated or removed sample, τ_(m) is the total time dilationexperienced up to t_(m), m is the number of sampling intervals removedor inserted to counteract time dilation, Δt is the sampling interval(time between successive ADC samples) and λ is the wavelength at thecentre frequency of the signal sampling band before any frequencydownconversion at a receiver 18. The parameter b′₁ is an estimate of thetime invariant component of DFO b₁ with any systematic error (due tosatellite turnaround frequency change) corrected using an estimatederived from the reference signal. Alternatively b′₁ may be anotherparameter which is varied in steps in maximising the modulus of the CCF.It is also possible to exploit b₂ in the above expression as anindependent parameter varied in steps, thereby dissociating timedilation and time variation of DFO.

[0092] If the effects of satellite motion are not too severe, it ispossible to omit the quadratic second term in Equation (14) in t_(m) ²,and apply the resulting linear correction.

[0093] Samples preceded by an omitted or inserted sample are multipliedby a phase correction factor P given by one of the following:

P=e ^(2πifΔt), or P=e ^(−2πifΔt)   (15)

[0094] The exponent in the expression for P is positive or negativeaccording to whether the preceding sample is omitted or inserted; f isthe centre frequency of the samples' signal after downconversion forinput to an analogue to digital converter, and Δt is the intervalbetween successive ADC samples. The samples' signal band is set by ananti-aliasing low-pass filter.

[0095] Before time dilation compensation is undertaken, the modulus ofthe CCF is maximised to obtain approximate values of DFO and theparameters τ₀, b₁ and b₂ as previously discussed. Time dilationcompensation is then implemented by inserting or excising samples asappropriate and adjusting phase in each of the four data sets z₁, z₂,p₁, p₂ to produce compensated equivalents. The CCF is then maximisedagain using the compensated data sets so formed to provide better valuesof DTO, DFO and if required DFRO and DDROA.

[0096] The process of compensation for change in both DFO motion andtime dilation is summarised as follows:

[0097] (a) obtaining initial values of DFO and DTO using compensationfor varying DFO only; (b) time dilation compensation of one data set bydelaying/advancing (translating) parts of it relative to an associateddata set and applying phase corrections; this involvesreplicating/deleting samples as appropriate; and

[0098] (c) multiplying the sample sets (in the sense of an innerproduct) to obtain a measure of the signal energy and to vary parametervalues to maximise the modulus of the CCF given in general terms by:

CCF=[z ₁ *]T ⁻¹ D ⁻¹ [z ₂]  (16)

[0099] It is the CCF modulus which is maximised because noise fromreceiving equipment is included, which makes signal energy a complexquantity (in the mathematical sense). The T⁻¹ and D⁻¹ notation is usedto indicate translation compensation and frequency distortioncompensation respectively, the superscripts denoting reversal of effectsintroduced by system geometry and kinematics. The brackets [ ]indicate avector comprising the signal samples. Translation of a relay satelliteleads to time dilation because it alters the signal path and thereforealso the associated time of flight.

[0100] In developing more general expressions for the DTO and DFO itmust be borne in mind that the optimisation of the resulting CCF becomesmore processor intensive. The approximation described herein using τ₀,b₁ and b₂ is well within the capabilities of current parallelprocessors.

[0101] It is expedient to represent the distortion compensation in termsof a DFO compensation and a separate DTO compensation because a rapidchange of phase is removed by this approach, leaving the DTO dilationcompensation acting on one data stream without the need forinterpolation between signal samples.

[0102] Referring now to FIGS. 8 to 11, the effects of counteracting DFOvariation and time dilation in accordance with the invention areillustrated graphically. These drawings are plots of DTO against DFOderived using simulated data without significant noise. In FIG. 8, DFOvaries with DTO over a region 80 which includes error limits. FIG. 9illustrates the effect of frequency compensation alone, ie counteractingchange in DFO. A region 90 is produced which is parallel to the DFO axisindicating DFO is now constant. FIG. 10 is the equivalent of FIG. 9 fortime or DTO compensation alone, and shows a region 100 of constant DTO.The combination of both these procedures is shown in FIG. 11, whichshows a single peak of constant DFO and DTO. Compensation for DTOvariation was carried out using the linear term in Equation (14) only.

[0103] FIGS. 12 to 14 illustrate graphically the effect of introducingnoise. Except as regards noise they were produced using the same data asFIGS. 8 to 11, and they would show precisely the same features absentthe noise. It can be seen that DTO and DFO are not visible in FIGS. 12and 13, either without compensation (FIG. 12) or with compensation foronly one of frequency distortion and time distortion (FIG. 13), whichgive like results. FIG. 14 has a single peak 140 easily discernableabove surrounding noise such as 142; it demonstrates that compensationfor both frequency and time distortion in accordance with this inventionleads to DTO and DFO values being obtainable despite the presence ofnoise.

1. A method of locating the source of an unknown signal received by aplurality of signal relays, the method including the steps of: (a)arranging for a plurality of receivers (18) to receive replicas of theunknown signal via respective signal relays (14, 16); and (b) subjectingthe replicas to correlation processing; characterised in thatcorrelation processing in step (b) is performed with a complexcorrelation function (CCF) at least partly compensated for change in thereplicas' Differential Frequency Offset (DFO) with time due to relaymotion relative to the source and receivers.
 2. A method according toclaim 1 characterised in that correlation processing in step (b)includes introducing a trial time offset between the replicas andevaluating their correlation and iterating this to obtain a correlationmaximum and derive at least one of the replicas' Differential TimeOffset (DTO) and DFO;
 3. A method according to claim 2 characterised inthat correlation processing in step (b) is carried out with a CCFcontaining an exponent of a function of time having a first term whichis a constant DFO value and a second term which is a product of time anda constant value for rate of change of DFO with time, ie a constantdifferential frequency rate offset (DFRO) value, and wherein step(b))also includes introducing a trial value corresponding to DFRO priorto evaluating correlation and iterating one type of trial value for eachof the other type and determining a DFRO value appropriate to obtain acorrelation maximum.
 4. A method according to claim 1, 2 or 3characterised in that correlation processing in step (b) is performedwith data sets adapted in phase and subject to data replication orremoval to counteract time dilation.
 5. A method according to claim 4characterised in that samples in one of two data sets to be correlatedare adapted in phase by samples following a removed sample beingmultiplied by a phase factor or by samples including and following areplicated sample being multiplied by such a factor, the phase factorbeing either or e^(2πifΔt) or e^(−2πifΔt) according to which of the twodata sets is adapted and to the sign of DTO change with time, andwherein f is a signal band centre frequency of the samples afterdownconversion for sampling and Δt is an interval between successivesamples.
 6. A method according to claim 4 or 5 characterised in thatsamples are selected for removal or replication as the case may be onthe basis of the correlation function's time dimension spread over aninterval in which the samples were obtained.
 7. A method according toclaim 4 or 5 characterised in that samples are selected for removal orreplication as the case may be on the basis of at least a linearapproximation to the following quadratic equation for t_(m):$\tau_{m} = {{{{- \frac{\lambda}{c}}b_{1}^{\prime}t_{m}} - {\frac{\lambda}{c}b_{2}\tau_{m}^{2}}} = {{m\Delta}\quad t}}$

where t_(m) is the time to the mth replicated or removed sample, m isthe number of replicated or removed samples up to time t_(m), τ_(m) isthe total time dilation up to time t_(m), Δt is the time intervalbetween successive samples and λ is the free space wavelength at thecentre frequency of the sampling band, b′₁ is an estimate of DFO atinitiation of sampling corrected for any systematic error in frequency,and b₂ is half of the rate of change of DFO with time.
 8. A methodaccording to any one of claims 1 to 7 characterised in that thecorrelation function is a complex correlation function referred to as aCCF given by:CCF ≡ A(τ₀, b₁, b₂) = ∫₀^(T)z₁^(*)(t)z(t + τ₀)^(−2π  i(b₁ + 2b₂t)t)  t

where z₁ and z₂ are data sets representing two signals being correlated,the asterisk denotes the complex conjugate of z₁, T is the time overwhich samples are taken, τ₀ is a time invariant component of DTO, and b₁and b₂ are constants in an expression for change in DFO with time t asfollows:—DFO≡ν=b₁+2b₂t, and wherein the CCF is maximised inter alia byinserting therein trial values of b₂ and versions of z₁ and z₂ adaptedto contain added or removed sample values to counteract time dilation.9. A method according to claim 8 characterised in that after eachiteration of step (b) the product of z₁*(t)z₂(t+τ₀)exp(−2πib₂t²) isFourier transformed to the frequency domain, in order to provide for DFOb₁ to be obtainable as the frequency offset at which the CCF maximumoccurs.
 10. A location device for locating the source (10) of an unknownsignal received by a plurality of signal relays (14, 16), and including:(a) a plurality of receivers (18) for receiving replicas of the unknownsignal via respective signal relays (14, 16); (b) a correlationprocessor (30) for correlating the replicas and obtaining a correlationmaximum indicating at least one of their DTO and DFO; characterised inthat the correlation processor (30) is arranged to perform correlationwith a complex correlation function (CCF) at least partly compensatedfor change in DFO with time due to relay motion relative to the sourceand receivers.
 11. A device according to claim 10 characterised in thatthe correlation processor (30) is arranged to introduce a trial timeoffset between the replicas, evaluate their correlation and iterate thisprocedure to obtain a correlation maximum.
 12. A device according toclaim 11 characterised in that the correlation processor (30) isarranged to employ a CCF which contains an exponent of a function oftime having a first term which is a constant DFO value and a second termwhich is a product of time and a constant value for rate of change ofDFO with time, ie a constant differential frequency rate offset (DFRO)value, and wherein the correlation processor (30) is also arranged tointroduce trial values of DFRO and to iteratively evaluate correlationfor pairs of trial values of DTO and DFRO to determine a DFRO valueenabling a correlation maximum to be obtained.
 13. A device according toclaim 11 or 12 characterised in that the correlation processor (30) isarranged to perform correlation with data sets adapted in phase andsubject to data replication or removal to counteract time dilation. 14.A device according to claim 13 characterised in that the correlationprocessor (30) is arranged to adapt in phase one of two sets of datasamples to be correlated by multiplying samples following a removedsample by a phase factor or multiplying a replicated sample and samplesfollowing it by a such a factor, the phase factor being eithere^(2πifΔt) or e^(−2πifΔt) according to which of the two data sets isadapted and to the sign of DTO change with time, and wherein f is asignal band centre frequency of the samples after downconversion forsampling and Δt is an interval between successive samples.
 15. A deviceaccording to claim 13 or 14 characterised in that the correlationprocessor (30) is arranged to select samples for removal and replicationon the basis of the correlation function's time dimension spread over aninterval in which the samples were obtained.
 16. A device according toclaim 13 or 14 characterised in that the correlation processor (30) isarranged to select samples for removal and replication on the basis ofat least the linear term of the following quadratic equation in t_(m) ²:$\tau_{m} = {{{{- \frac{\lambda}{c}}b_{1}^{\prime}t_{m}} - {\frac{\lambda}{c}b_{2}\tau_{m}^{2}}} = {{m\Delta}\quad t}}$

where t_(m) is the time from initiation of sampling to the mthreplicated or removed sample, m is the number of replicated or removedsamples up to time t_(m), τ_(m) is the total time dilation up to timet_(m), Δt is the time interval between successive samples, λ is the freespace wavelength at the centre frequency of the sampling band, b′₁ is anestimate of DFO at initiation of sampling corrected for any systematicerror in frequency, and b₂ is half of the rate of change of DFO withtime.
 17. A device according to claim 10 characterised in that the CCFis given by:CCF ≡ A(τ₀, b₁, b₂) = ∫₀^(T)z₁^(*)(t)z(t + τ₀)^(−2π  i(b₁ + 2b₂t)t)  t

where z₁ and z₂ are data sets representing two signals being correlated,the asterisk denotes the complex conjugate of z₁, T is the time overwhich samples are taken, τ₀ is a time invariant component of DTO, and b₁and b₂ are constants in an expression for DFO for change in DFO withtime t as follows:—DFO≡ν=b₁+2b₂t, and wherein the correlation processor(30) is arranged to maximise the CCF by evaluating it with trial valuesof b₂ in addition to those of τ₀, together with versions of z₁ and z₂adapted to contain replicated or removed sample values to counteracttime dilation.
 18. A device according to claim 17 characterised in thatthe correlation processor (30) is arranged to Fourier transformsuccessive products z₁*(t)z₂(t+τ₀)exp(−2πib₂t²) to the frequency domainin order to obtain DFO b₁ as a frequency offset at which a CCF maximumoccurs.